//
// Created by Ivor_ on 2022/1/22.
//

#include "matrixCalculate.h"

double* dMatLVec(int* MatDims, double** mat, double* vec) {
    double* retVec = dNewV(0).create(MatDims[1]);
    for (int i = 0; i < MatDims[1]; ++i) {
        for (int j = 0; j < MatDims[0]; ++j) {
            retVec[i] += (vec[j]*mat[j][i]);
        }
    }
    return retVec;
}

double* dMatRVec(int* MatDims, double** mat, double* vec) {
    double* retVec = dNewV(0).create(MatDims[0]);
    for (int i = 0; i < MatDims[0]; ++i) {
        for (int j = 0; j < MatDims[1]; ++j) {
            retVec[i] += (mat[i][j]*vec[j]);
        }
    }
    return retVec;
}

_Complex double* cMatLVec(int* MatDims, _Complex double** mat, _Complex double* vec) {
    _Complex double* retVec = cNewV(0).create(MatDims[1]);
    for (int i = 0; i < MatDims[1]; ++i) {
        for (int j = 0; j < MatDims[0]; ++j) {
            retVec[i] = retVec[i] + (vec[j]*mat[j][i]);
        }
    }
    return retVec;
}

_Complex double* cMatRVec(int* MatDims, _Complex double** mat, _Complex double* vec) {
    _Complex double* retVec = cNewV(0).create(MatDims[0]);
    for (int i = 0; i < MatDims[0]; ++i) {
        for (int j = 0; j < MatDims[1]; ++j) {
            retVec[i] = retVec[i] + (mat[i][j]*vec[j]);
        }
    }
    return retVec;
}

double dMixVal(int* MatDims, double** mat, double* vecL, double* vecR) {
    return dvf().innerProd(MatDims[1], dMatLVec(MatDims, mat, vecL), vecR);
}

_Complex double cMixVal(int* MatDims, _Complex double** mat, _Complex double* vecL, _Complex double* vecR) {
    return cvf().innerProd(MatDims[1], cMatLVec(MatDims, mat, vecL), vecR);
}

double dQuadratic(int* MatDims, double** mat, double* vec) {
    return dMixVal(MatDims, mat, vec, vec);
}

_Complex double cQuadratic(int* MatDims, _Complex double** mat, _Complex double* vec) {
    return cMixVal(MatDims, mat, vec, vec);
}

double dRayleighQuo(int* MatDims, double** mat, double* vec) {
    return (dQuadratic(MatDims, mat, vec) / dvf().innerProd(MatDims[0], vec, vec));
}

_Complex double cRayleighQuo(int* MatDims, _Complex double** mat, _Complex double* vec) {
    return cQuadratic(MatDims, mat, vec)/cvf().innerProd(MatDims[0], vec, vec);
}

double dTrace(int dim, double** matrix) {
    double retVal = 0.;
    for (int i = 0; i < dim; ++i) {
        retVal += matrix[i][i];
    }
    return retVal;
}

_Complex double cTrace(int dim, _Complex double** matrix) {
    _Complex double retVal = 0. + 0.i;
    for (int i = 0; i < dim; ++i) {
        retVal += matrix[i][i];
    }
    return retVal;
}

double** * dNaiveMatMultiple(int* dims, double** matL, double** matR) {
    double** * retMat = (double***)malloc(sizeof(double**));
    *retMat = dNewM(SHAPE_RECTANGLE, 0).create_rm(simDim(2, dims[0], dims[3]));
    for (int i = 0; i < dims[0]; ++i) {
        for (int j = 0; j < dims[3]; ++j) {
            for (int k = 0; k < dims[1]; ++k) {
                (*retMat)[i][j] += (matL[i][k]*matR[k][j]);
            }
        }
    }
    return retMat;
}

_Complex double** * cNaiveMatMultiple(int* dims, _Complex double** matL, _Complex double** matR) {
    _Complex double** * retMat = (_Complex double***)malloc(sizeof(_Complex double**));
    *retMat = cNewM(SHAPE_RECTANGLE, 0).create_rm(simDim(2, dims[0], dims[3]));
    for (int i = 0; i < dims[0]; ++i) {
        for (int j = 0; j < dims[3]; ++j) {
            for (int k = 0; k < dims[1]; ++k) {
                (*retMat)[i][j] += (matL[i][k]*matR[k][j]);
            }
        }
    }
    return retMat;
}

double** * dNaiveMatHadamardMul(int* dims, double** mat1, double** mat2) {
    double** * retMat = (double***)malloc(sizeof(double**));
    *retMat = dNewM(SHAPE_RECTANGLE, 0).create_rm(simDim(2, dims[0], dims[2]));
    for (int i = 0; i < dims[0]; ++i) {
        for (int j = 0; j < dims[1]; ++j) {
            (*retMat)[i][j] = mat1[i][j]*mat2[i][j];
        }
    }
    return retMat;
}

_Complex double** * cNaiveMatHadamardMul(int* dims, _Complex double** mat1, _Complex double** mat2) {
    _Complex double** * retMat = (_Complex double***)malloc(sizeof(_Complex double**));
    *retMat = cNewM(SHAPE_RECTANGLE, 0).create_rm(simDim(2, dims[0], dims[2]));
    for (int i = 0; i < dims[0]; ++i) {
        for (int j = 0; j < dims[1]; ++j) {
            (*retMat)[i][j] = mat1[i][j]*mat2[i][j];
        }
    }
    return retMat;
}

double** * dNaiveKroneckerMul(int* dims, double** matL, double** matR) {
    double** * retMat = (double***)malloc(sizeof(double**));
    *retMat = dNewM(SHAPE_RECTANGLE, 0).create_rm(simDim(2, dims[0]*dims[2], dims[1]*dims[3]));
    for (int i = 0; i < dims[0]*dims[2]; ++i) {
        for (int j = 0; j < dims[1]*dims[3]; ++j) {
            (*retMat)[i][j] = matL[i/dims[2]][j/dims[3]] * matR[i%dims[2]][j%dims[3]];
        }
    }
    return retMat;
}

_Complex double** * cNaiveKroneckerMul(int* dims, _Complex double** matL, _Complex double** matR) {
    _Complex double** * retMat = (_Complex double***)malloc(sizeof(_Complex double**));
    *retMat = cNewM(SHAPE_RECTANGLE, 0).create_rm(simDim(2, dims[0]*dims[2], dims[1]*dims[3]));
    for (int i = 0; i < dims[0]*dims[2]; ++i) {
        for (int j = 0; j < dims[1]*dims[3]; ++j) {
            (*retMat)[i][j] = matL[i/dims[2]][j/dims[3]]*matR[i%dims[2]][j%dims[3]];
        }
    }
    return retMat;
}

double** * dNaiveMatInnMul(int* dims, double** matL, double** matR) {
    return dNaiveMatMultiple(simDim(4, dims[1], dims[0], dims[2], dims[3]), dCvtM().trans(simDim(2, dims[0], dims[1]), matL), matR);
}

_Complex double** * cNaiveMatInnMul(int* dims, _Complex double** matL, _Complex double** matR) {
    return cNaiveMatMultiple(simDim(4, dims[1], dims[0], dims[2], dims[3]), cCvtM().dag(simDim(2, dims[0], dims[1]), matL), matR);
}

double** * dNaiveMatOutMul(int* dims, double** matL, double** matR) {
    return dNaiveMatMultiple(simDim(4, dims[0], dims[1], dims[3], dims[2]), matL, dCvtM().trans(simDim(2, dims[2], dims[3]), matR));
}

_Complex double** * cNaiveMatOutMul(int* dims, _Complex double** matL, _Complex double** matR) {
    return cNaiveMatMultiple(simDim(4, dims[0], dims[1], dims[3], dims[2]), matL, cCvtM().dag(simDim(2, dims[2], dims[3]), matR));
}

double** * dUltraMatMultiple(int* dims, int count, ...) {
    va_list vl;
    va_start(vl, count);
    double** * retMat = (double***)malloc(sizeof(double**));
    *retMat = dNewM(SHAPE_RECTANGLE, 0).create_rm(simDim(2, dims[0], dims[1]));
    double** tmp = va_arg(vl, double**);
    for (int i = 0; i < dims[0]; ++i) {
        for (int j = 0; j < dims[1]; ++j) {
            (*retMat)[i][j] = tmp[i][j];
        }
    }
    for (int i = 0; i < count-1; ++i) {
        double** * cmt = dNaiveMatMultiple(simDim(4, dims[0], dims[2*i+1], dims[2*i+2], dims[2*i+3]), *retMat, va_arg(vl, double**));
        for (int j = 0; j < dims[0]; ++j) {
            for (int k = 0; k < dims[2*i+3]; ++k) {
                (*retMat)[j][k] = (*cmt)[j][k];
            }
        }
        freeDMat(dims[0], *cmt);
    }
    va_end(vl);
    return retMat;
}

_Complex double** * cUltraMatMultiple(int* dims, int count, ...) {
    va_list vl;
    va_start(vl, count);
    _Complex double** * retMat = (_Complex double***)malloc(sizeof(_Complex double**));
    *retMat = cNewM(SHAPE_RECTANGLE, 0).create_rm(simDim(2, dims[0], dims[1]));
    _Complex double** tmp = va_arg(vl, _Complex double**);
    for (int i = 0; i < dims[0]; ++i) {
        for (int j = 0; j < dims[1]; ++j) {
            (*retMat)[i][j] = tmp[i][j];
        }
    }
    for (int i = 0; i < count-1; ++i) {
        _Complex double** * cmt = cNaiveMatMultiple(simDim(4, dims[0], dims[2*i+1], dims[2*i+2], dims[2*i+3]), *retMat, va_arg(vl, _Complex double**));
        for (int j = 0; j < dims[0]; ++j) {
            for (int k = 0; k < dims[2*i+3]; ++k) {
                (*retMat)[j][k] = (*cmt)[j][k];
            }
        }
        freeCMat(dims[0], *cmt);
    }
    va_end(vl);
    return retMat;
}

double** * dUltraMatHadamardMul(int* dims, int count, ...) {
    va_list vl;
    va_start(vl, count);
    double** * retMat = (double***)malloc(sizeof(double**));
    *retMat = dNewM(SHAPE_RECTANGLE, 0).create_rm(simDim(2, dims[0], dims[1]));
    double** tmp = va_arg(vl, double**);
    for (int i = 0; i < dims[0]; ++i) {
        for (int j = 0; j < dims[1]; ++j) {
            (*retMat)[i][j] = tmp[i][j];
        }
    }
    for (int i = 0; i < count-1; ++i) {
        double** * cmt = dNaiveMatHadamardMul(simDim(4, dims[0], dims[1], dims[2*i+2], dims[2*i+3]), *retMat, va_arg(vl, double**));
        for (int j = 0; j < dims[0]; ++j) {
            for (int k = 0; k < dims[1]; ++k) {
                (*retMat)[j][k] = (*cmt)[j][k];
            }
        }
        freeDMat(dims[0], *cmt);
    }
    va_end(vl);
    return retMat;
}

_Complex double** * cUltraMatHadamardMul(int* dims, int count, ...) {
    va_list vl;
    va_start(vl, count);
    _Complex double** * retMat = (_Complex double***)malloc(sizeof(_Complex double**));
    *retMat = cNewM(SHAPE_RECTANGLE, 0).create_rm(simDim(2, dims[0], dims[1]));
    _Complex double** tmp = va_arg(vl, _Complex double**);
    for (int i = 0; i < dims[0]; ++i) {
        for (int j = 0; j < dims[1]; ++j) {
            (*retMat)[i][j] = tmp[i][j];
        }
    }
    for (int i = 0; i < count-1; ++i) {
        _Complex double** * cmt = cNaiveMatHadamardMul(simDim(4, dims[0], dims[1], dims[2*i+2], dims[2*i+3]), *retMat, va_arg(vl, _Complex double**));
        for (int j = 0; j < dims[0]; ++j) {
            for (int k = 0; k < dims[1]; ++k) {
                (*retMat)[j][k] = (*cmt)[j][k];
            }
        }
        freeCMat(dims[0], *cmt);
    }
    va_end(vl);
    return retMat;
}

double** * dUltraKroneckerMul(int* dims, int count, ...) {
    va_list vl;
    int mc = 1, mr = 1;
    int cc = 1, cr = 1;
    for (int i = 0; i < count; ++i) {
        mc *= dims[i*2];
        mr *= dims[i*2+1];
    }
    va_start(vl, count);
    double** * retMat = (double***)malloc(sizeof(double**));
    *retMat = dNewM(SHAPE_RECTANGLE, 0).create_rm(simDim(2, mc, mr));
    double** tmp = va_arg(vl, double**);
    for (int i = 0; i < dims[0]; ++i) {
        for (int j = 0; j < dims[1]; ++j) {
            (*retMat)[i][j] = tmp[i][j];
        }
    }
    for (int i = 0; i < count-1; ++i) {
        cc *= dims[2*i]; cr *= dims[2*i+1];
        double** * cmt = dNaiveKroneckerMul(simDim(4, cc, cr, dims[2*i+2], dims[2*i+3]), *retMat, va_arg(vl, double**));
        for (int j = 0; j < cc*dims[2*i+2]; ++j) {
            for (int k = 0; k < cr*dims[2*i+3]; ++k) {
                (*retMat)[j][k] = (*cmt)[j][k];
            }
        }
        freeDMat(cc*dims[2*i+2], *cmt);
    }
    va_end(vl);
    return retMat;
}

_Complex double** * cUltraKroneckerMul(int* dims, int count, ...) {
    va_list vl;
    int mc = 1, mr = 1;
    int cc = 1, cr = 1;
    for (int i = 0; i < count; ++i) {
        mc *= dims[i*2];
        mr *= dims[i*2+1];
    }
    va_start(vl, count);
    _Complex double** * retMat = (_Complex double***)malloc(sizeof(_Complex double**));
    *retMat = cNewM(SHAPE_RECTANGLE, 0).create_rm(simDim(2, mc, mr));
    _Complex double** tmp = va_arg(vl, _Complex double**);
    for (int i = 0; i < dims[0]; ++i) {
        for (int j = 0; j < dims[1]; ++j) {
            (*retMat)[i][j] = tmp[i][j];
        }
    }
    for (int i = 0; i < count-1; ++i) {
        cc *= dims[2*i]; cr *= dims[2*i+1];
        _Complex double** * cmt = cNaiveKroneckerMul(simDim(4, cc, cr, dims[2*i+2], dims[2*i+3]), *retMat, va_arg(vl, _Complex double**));
        for (int j = 0; j < cc*dims[2*i+2]; ++j) {
            for (int k = 0; k < cr*dims[2*i+3]; ++k) {
                (*retMat)[j][k] = (*cmt)[j][k];
            }
        }
        freeCMat(cc*dims[2*i+2], *cmt);
    }
    va_end(vl);
    return retMat;
}

double** * dUltraMatInnMul(int* dims, int count, ...) {
    va_list vl;
    va_start(vl, count);
    double** * retMat = (double***)malloc(sizeof(double**));
    *retMat = dNewM(SHAPE_RECTANGLE, 0).create_rm(simDim(2, dims[0], dims[1]));
    double** tmp = va_arg(vl, double**);
    for (int i = 0; i < dims[0]; ++i) {
        for (int j = 0; j < dims[1]; ++j) {
            *retMat[i][j] = tmp[i][j];
        }
    }
    for (int i = 0; i < count-1; ++i) {
        double** * cmt = dNaiveMatInnMul(simDim(4, dims[i==0 ? 0 : 2*i-1], dims[2*i+1], dims[2*i+2], dims[2*i+3]), *retMat, va_arg(vl, double**));
        for (int j = 0; j < dims[2*i+1]; ++j) {
            for (int k = 0; k < dims[2*i+3]; ++k) {
                (*retMat)[j][k] = (*cmt)[j][k];
            }
        }
        freeDMat(dims[2*i+1], *cmt);
    }
    va_end(vl);
    return retMat;
}

_Complex double** * cUltraMatInnMul(int* dims, int count, ...) {
    va_list vl;
    va_start(vl, count);
    _Complex double** * retMat = (_Complex double***)malloc(sizeof(_Complex double**));
    (*retMat) = cNewM(SHAPE_RECTANGLE, 0).create_rm(simDim(2, dims[0], dims[1]));
    _Complex double** tmp = va_arg(vl, _Complex double**);
    for (int i = 0; i < dims[0]; ++i) {
        for (int j = 0; j < dims[1]; ++j) {
            (*retMat)[i][j] = tmp[i][j];
        }
    }
    for (int i = 0; i < count-1; ++i) {
        _Complex double** * cmt = cNaiveMatInnMul(simDim(4, dims[i==0 ? 0 : 2*i-1], dims[2*i+1], dims[2*i+2], dims[2*i+3]), *retMat, va_arg(vl, _Complex double**));
        for (int j = 0; j < dims[2*i+1]; ++j) {
            for (int k = 0; k < dims[2*i+3]; ++k) {
                (*retMat)[j][k] = (*cmt)[j][k];
            }
        }
        freeCMat(dims[2*i+1], *cmt);
    }
    va_end(vl);
    return retMat;
}
double** * dUltraMatOutMul(int* dims, int count, ...) {
    va_list vl;
    va_start(vl, count);
    double** * retMat = (double***)malloc(sizeof(double**));
    *retMat = dNewM(SHAPE_RECTANGLE, 0).create_rm(simDim(2, dims[0], dims[1]));
    double** tmp = va_arg(vl, double**);
    for (int i = 0; i < dims[0]; ++i) {
        for (int j = 0; j < dims[1]; ++j) {
            (*retMat)[i][j] = tmp[i][j];
        }
    }
    for (int i = 0; i < count-1; ++i) {
        double** * cmt = dNaiveMatOutMul(simDim(4, dims[0], dims[i==0 ? 1 : 2*i], dims[2*i+2], dims[2*i+3]), *retMat, va_arg(vl, double**));
        for (int j = 0; j < dims[0]; ++j) {
            for (int k = 0; k < dims[2*i+2]; ++k) {
                (*retMat)[j][k] = (*cmt)[j][k];
            }
        }
        freeDMat(dims[0], *cmt);
    }
    va_end(vl);
    return retMat;
}

_Complex double** * cUltraMatOutMul(int* dims, int count, ...) {
    va_list vl;
    va_start(vl, count);
    _Complex double** * retMat = (_Complex double***)malloc(sizeof(_Complex double**));
    *retMat = cNewM(SHAPE_RECTANGLE, 0).create_rm(simDim(2, dims[0], dims[1]));
    _Complex double** tmp = va_arg(vl, _Complex double**);
    for (int i = 0; i < dims[0]; ++i) {
        for (int j = 0; j < dims[1]; ++j) {
            (*retMat)[i][j] = tmp[i][j];
        }
    }
    for (int i = 0; i < count-1; ++i) {
        _Complex double** * cmt = cNaiveMatOutMul(simDim(4, dims[0], dims[i==0 ? 1 : 2*i], dims[2*i+2], dims[2*i+3]), *retMat, va_arg(vl, _Complex double**));
        for (int j = 0; j < dims[0]; ++j) {
            for (int k = 0; k < dims[2*i+2]; ++k) {
                (*retMat)[j][k] = (*cmt)[j][k];
            }
        }
        freeCMat(dims[0], *cmt);
    }
    va_end(vl);
    return retMat;
}

struct dMatrixNaiveMultiple dMulN(void) {
    struct dMatrixNaiveMultiple dnm;
    dnm.multiple = dNaiveMatMultiple;
    dnm.hadamard = dNaiveMatHadamardMul;
    dnm.kronecker = dNaiveKroneckerMul;
    dnm.inner = dNaiveMatInnMul;
    dnm.outer = dNaiveMatOutMul;
    return dnm;
}

struct cMatrixNaiveMultiple cMulN(void) {
    struct cMatrixNaiveMultiple cnm;
    cnm.multiple = cNaiveMatMultiple;
    cnm.hadamard = cNaiveMatHadamardMul;
    cnm.kronecker = cNaiveKroneckerMul;
    cnm.inner = cNaiveMatInnMul;
    cnm.outer = cNaiveMatOutMul;
    return cnm;
}

struct dMatrixUltraMultiple dMulX(void) {
    struct dMatrixUltraMultiple dxm;
    dxm.multiple = dUltraMatMultiple;
    dxm.hadamard = dUltraMatHadamardMul;
    dxm.kronecker = dUltraKroneckerMul;
    dxm.inner = dUltraMatInnMul;
    dxm.outer = dUltraMatOutMul;
    return dxm;
}

struct cMatrixUltraMultiple cMulX(void) {
    struct cMatrixUltraMultiple cxm;
    cxm.multiple = cUltraMatMultiple;
    cxm.hadamard = cUltraMatHadamardMul;
    cxm.kronecker = cUltraKroneckerMul;
    cxm.inner = cUltraMatInnMul;
    cxm.outer = cUltraMatOutMul;
    return cxm;
}

void freeDMat(int lines, double** mat) {
    for (int i = 0; i < lines; ++i) {
        free(mat[i]);
    }
    free(mat);
}

void freeCMat(int lines, _Complex double** mat) {
    for (int i = 0; i < lines; ++i) {
        free(mat[i]);
    }
    free(mat);
}

void freeDVec(double* vec) {
    free(vec);
}

void freeCVec(_Complex double* vec) {
    free(vec);
}

void dSchmidt(int dim, double** matrix)  {
    double temp = 0.;
    for (int i = 0; i < dim; ++i) {
        temp += matrix[i][0]*matrix[i][0];
    }
    for (int i = 0; i < dim; ++i) {
        matrix[i][0] /= sqrt(temp);
    }
    for (int i = 1; i < dim; ++i) {
        double* tmpVec = (double*)calloc(dim, sizeof(double));
        for (int k = 0; k < dim; ++k) {
            tmpVec[k] = 0.;
        }
        for (int j = 0; j < i; ++j) {
            temp = 0.;
            for (int k = 0; k < dim; ++k) {
                temp += matrix[k][j]*matrix[k][i];
            }
            for (int k = 0; k < dim; ++k) {
                tmpVec[k] += (temp*matrix[k][j]);
            }
        }
        for (int k = 0; k < dim; ++k) {
            matrix[k][i] -= tmpVec[k];
        }
        freeDVec(tmpVec);
        temp = 0.;
        for (int j = 0; j < dim; ++j) {
            temp += matrix[j][i]*matrix[j][i];
        }
        for (int k = 0; k < dim; ++k) {
            matrix[k][i] /= sqrt(temp);
        }
    }
}

void cSchmidt(int dim, _Complex double** matrix)  {
    _Complex double temp = 0. + 0.i;
    for (int i = 0; i < dim; ++i) {
        temp += conj(matrix[i][0])*matrix[i][0];
    }
    for (int i = 0; i < dim; ++i) {
        matrix[i][0] /= sqrt(creal(temp));
    }
    for (int i = 1; i < dim; ++i) {
        _Complex double* tmpVec = (_Complex double*)calloc(dim, sizeof(_Complex double));
        for (int k = 0; k < dim; ++k) {
            tmpVec[k] = 0. + 0.i;
        }
        for (int j = 0; j < i; ++j) {
            temp = 0. + 0.i;
            for (int k = 0; k < dim; ++k) {
                temp += conj(matrix[k][j])*matrix[k][i];
            }
            for (int k = 0; k < dim; ++k) {
                tmpVec[k] += (temp*matrix[k][j]);
            }
        }
        for (int k = 0; k < dim; ++k) {
            matrix[k][i] -= tmpVec[k];
        }
        freeCVec(tmpVec);
        temp = 0.;
        for (int j = 0; j < dim; ++j) {
            temp += conj(matrix[j][i])*matrix[j][i];
        }
        for (int k = 0; k < dim; ++k) {
            matrix[k][i] /= sqrt(temp);
        }
    }
}

double dDeterminant(int dim, double** matrix) {
    double** cpyMat = (double**)calloc(dim, sizeof(double*));
    for (int i = 0; i < dim; ++i) {
        cpyMat[i] = (double*)calloc(dim, sizeof(double));
        for (int j = 0; j < dim; ++j) {
            cpyMat[i][j] = matrix[i][j];
        }
    }
    int main = -1, signs = 0;
    double dmt = 1.;
    for (int j = 0; j < dim; j++) {
        int flag = 0;
        int mlo = j;
        for (int i = j; i < dim; i++) {
            if (cpyMat[i][j] != 0) {
                flag = 1;
                mlo = i;
                main++;
                break;
            }
        }
        if (flag == 1) {
            if (mlo != main) {
                dSwtM('l').exc(cpyMat, simDim(2, dim, dim), simDim(2, mlo, main));
                signs++;
            }
            for (int k = 0; k < dim; k++) {
                if (k != main) {
                    dSwtM('l').dbl(cpyMat, simDim(2, dim, dim), simDim(2, main, k), -cpyMat[k][j]/cpyMat[main][j]);
                }
            }
        }
    }
    for (int i = 0; i < dim; i++) {
        dmt *= cpyMat[i][i];
    }
    freeDMat(dim, cpyMat);
    dmt *= pow(-1.0, (double)signs);
    return dmt;
}

_Complex double cDeterminant(int dim, _Complex double** matrix) {
    _Complex double** cpyMat = (_Complex double**)calloc(dim, sizeof(_Complex double*));
    for (int i = 0; i < dim; ++i) {
        cpyMat[i] = (_Complex double*)calloc(dim, sizeof(_Complex double));
        for (int j = 0; j < dim; ++j) {
            cpyMat[i][j] = matrix[i][j];
        }
    }
    int main = -1, signs = 0;
    _Complex double dmt = 1. + 0.i;
    for (int j = 0; j < dim; j++) {
        int flag = 0;
        int mlo = j;
        for (int i = j; i < dim; i++) {
            if (cpyMat[i][j] == 0. + 0.i) {
                flag = 1;
                mlo = i;
                main++;
                break;
            }
        }
        if (flag == 1) {
            if (mlo != main) {
                cSwtM('l').exc(cpyMat, simDim(2, dim, dim), simDim(2, mlo, main));
                signs++;
            }
            for (int k = 0; k < dim; k++) {
                if (k != main) {
                    cSwtM('l').dbl(cpyMat, simDim(2, dim, dim), simDim(2, main, k), -(cpyMat[k][j])/cpyMat[main][j]);
                }
            }
        }
    }
    for (int i = 0; i < dim; i++) {
        dmt = dmt*cpyMat[i][i];
    }
    freeCMat(dim, cpyMat);
    dmt *= (pow(-1.0, (double)signs) + 0.i);
    return dmt;
}

double** dInverse(int dim, double** matrix) {
    int main = -1;
    double** cpyMat = (double**)calloc(dim, sizeof(double*));
    for (int i = 0; i < dim; ++i) {
        cpyMat[i] = (double*)calloc(dim, sizeof(double));
        for (int j = 0; j < dim; ++j) {
            cpyMat[i][j] = matrix[i][j];
        }
    }
    double** invMat = dNewM(SHAPE_SQUARE, 'U').create_sm(dim);
    for (int j = 0; j < dim; j++) {
        int flag = 0, mlo = j;
        for (int i = j; i < dim; i++) {
            if (cpyMat[i][j] != 0) {
                flag = 1;
                mlo = i;
                main++;
                break;
            }
        }
        if (flag == 1) {
            if (mlo != main) {
                dSwtM('l').exc(cpyMat, simDim(2, dim, dim), simDim(2, mlo, main));
                dSwtM('l').exc(invMat, simDim(2, dim, dim), simDim(2, mlo, main));
            }
            for (int k = 0; k < dim; k++) {
                if (k != main) {
                    double value = -cpyMat[k][j]/cpyMat[main][j];
                    dSwtM('l').dbl(cpyMat, simDim(2, dim, dim), simDim(2, main, k), value);
                    dSwtM('l').dbl(invMat, simDim(2, dim, dim), simDim(2, main, k), value);
                }
            }
        }
    }
    for (int i = 0; i < dim; i++) {
        dSwtM('l').mul(invMat, simDim(2, dim, dim), i, 1./cpyMat[i][i]);
        dSwtM('l').mul(cpyMat, simDim(2, dim, dim), i, 1./cpyMat[i][i]);
    }
    freeDMat(dim, cpyMat);
    return invMat;
}

_Complex double** cInverse(int dim, _Complex double** matrix) {
    int main = -1;
    _Complex double** cpyMat = (_Complex double**)calloc(dim, sizeof(_Complex double*));
    for (int i = 0; i < dim; ++i) {
        cpyMat[i] = (_Complex double*)calloc(dim, sizeof(_Complex double));
        for (int j = 0; j < dim; ++j) {
            cpyMat[i][j] = matrix[i][j];
        }
    }
    _Complex double** invMat = cNewM(SHAPE_SQUARE, 'U').create_sm(dim);
    for (int j = 0; j < dim; j++) {
        int flag = 0, mlo = j;
        for (int i = j; i < dim; i++) {
            if (cpyMat[i][j] == 0. + 0.i) {
                flag = 1;
                mlo = i;
                main++;
                break;
            }
        }
        if (flag == 1) {
            if (mlo != main) {
                cSwtM('l').exc(cpyMat, simDim(2, dim, dim), simDim(2, mlo, main));
                cSwtM('l').exc(invMat, simDim(2, dim, dim), simDim(2, mlo, main));
            }
            for (int k = 0; k < dim; k++) {
                if (k != main) {
                    _Complex double value = -cpyMat[k][j]/cpyMat[main][j];
                    cSwtM('l').dbl(cpyMat, simDim(2, dim, dim), simDim(2, main, k), value);
                    cSwtM('l').dbl(invMat, simDim(2, dim, dim), simDim(2, main, k), value);
                }
            }
        }
    }
    for (int i = 0; i < dim; i++) {
        cSwtM('l').mul(invMat, simDim(2, dim, dim), i, cpow(cpyMat[i][i], -1. + 0.i));
        cSwtM('l').mul(cpyMat, simDim(2, dim, dim), i, cpow(cpyMat[i][i], -1. + 0.i));
    }
    freeCMat(dim, cpyMat);
    return invMat;
}

double* dDiagJacobi(int dim, double** matrix, double** eigenVectors, double eps) {
    double sum = 0.;
    for (int i = 0; i < dim; ++i) {
        for (int j = 0; j < dim; ++j) {
            if (i != j) {
                sum += sqrt(pow(matrix[i][j], 2.));
            }
            eigenVectors[i][j] = 1.*(i == j);
        }
    }
    sum /= (dim*dim);
    double** cpyMat = (double**)calloc(dim, sizeof(double*));
    for (int i = 0; i < dim; ++i) {
        cpyMat[i] = (double*)calloc(dim, sizeof(double));
        for (int j = 0; j < dim; ++j) {
            cpyMat[i][j] = matrix[i][j];
        }
    }
    while (sum > eps) {
        int* loc = (int*)calloc(2, sizeof(int));
        loc[0] = 0, loc[1] = 1;
        for (int i = 0; i < dim-1; ++i) {
            for (int j = i+1; j < dim; ++j) {
                if (fabs(cpyMat[i][j]) > fabs(cpyMat[loc[0]][loc[1]])) {
                    loc[0] = i;
                    loc[1] = j;
                }
            }
        }
        double temp = cpyMat[loc[0]][loc[0]] - cpyMat[loc[1]][loc[1]];
        temp /= (2. * cpyMat[loc[1]][loc[0]]);
        double tanT = (cpyMat[loc[0]][loc[0]] == cpyMat[loc[1]][loc[1]] ? 1. : -temp - sqrt(pow(temp, 2.) + 1.));
        double cosT = 1. / sqrt(pow(tanT, 2.) + 1.);
        double sinT = tanT / sqrt(pow(tanT, 2.) + 1.);
        double** iterMat = (double**)calloc(dim, sizeof(double*));
        for (int i = 0; i < dim; ++i) {
            iterMat[i] = (double*)calloc(dim, sizeof(double));
            for (int j = 0; j < dim; ++j) {
                iterMat[i][j] = 0.;
            }
        }
        for (int i = 0; i < dim; ++i) {
            for (int j = 0; j < dim; ++j) {
                if (i == loc[0] && j == loc[0])
                    iterMat[i][j] = cpyMat[loc[0]][loc[0]] * pow(cosT, 2.) + cpyMat[loc[1]][loc[1]] * pow(sinT, 2.) + cpyMat[loc[0]][loc[1]] * 2. * sinT * cosT;
                else if (i == loc[1] && j == loc[1])
                    iterMat[i][j] = cpyMat[loc[0]][loc[0]] * pow(sinT, 2.) + cpyMat[loc[1]][loc[1]] * pow(cosT, 2.) - cpyMat[loc[0]][loc[1]] * 2. * sinT * cosT;
                else if ((i == loc[0] && j == loc[1]) || (i == loc[1] && j == loc[0]))
                    iterMat[i][j] = cpyMat[loc[0]][loc[1]] * (pow(cosT, 2.) - pow(sinT, 2.)) + 0.5 * (cpyMat[loc[1]][loc[1]] - cpyMat[loc[0]][loc[0]]) * 2. * sinT * cosT;
                else if (i == loc[0])
                    iterMat[i][j] = cpyMat[j][loc[0]] * cosT + cpyMat[j][loc[1]] * sinT;
                else if (j == loc[0])
                    iterMat[i][j] = cpyMat[i][loc[0]] * cosT + cpyMat[i][loc[1]] * sinT;
                else if (i == loc[1])
                    iterMat[i][j] = -cpyMat[j][loc[0]] * sinT + cpyMat[j][loc[1]] * cosT;
                else if (j == loc[1])
                    iterMat[i][j] = -cpyMat[i][loc[0]] * sinT + cpyMat[i][loc[1]] * cosT;
                else
                    iterMat[i][j] = cpyMat[i][j];
            }
        }
        for (int i = 0; i < dim; ++i) {
            for (int j = 0; j < dim; ++j) {
                cpyMat[i][j] = iterMat[i][j];
            }
        }
        freeDMat(dim, iterMat);
        sum = 0.;
        for (int i = 0; i < dim; ++i) {
            for (int j = 0; j < dim; ++j) {
                if (i != j) {
                    sum += sqrt(pow(cpyMat[i][j], 2.));
                }
            }
        }
        sum /= (dim*dim);
        double** vt = dNewM(SHAPE_SQUARE, 'U').create_sm(dim);
        vt[loc[0]][loc[0]] = cosT;
        vt[loc[1]][loc[1]] = cosT;
        vt[loc[0]][loc[1]] = sinT;
        vt[loc[1]][loc[0]] = -sinT;
        free(loc);
        iterMat = (double**)calloc(dim, sizeof(double*));
        for (int k = 0; k < dim; ++k) {
            iterMat[k] = (double*)calloc(dim, sizeof(double));
            for (int l = 0; l < dim; ++l) {
                iterMat[k][l] = 0.;
                for (int m = 0; m < dim; ++m) {
                    iterMat[k][l] += vt[k][m]*eigenVectors[m][l];
                }
            }
        }
        for (int k = 0; k < dim; ++k) {
            for (int l = 0; l < dim; ++l) {
                eigenVectors[k][l] = iterMat[k][l];
            }
        }
        freeDMat(dim, iterMat);
    }
    double* eigenvalues = (double*)calloc(dim, sizeof(double));
    for (int i = 0; i < dim; ++i) {
        eigenvalues[i] = cpyMat[i][i];
    }
    freeDMat(dim, cpyMat);
    double** tempMat = (double**)calloc(dim, sizeof(double*));
    for (int k = 0; k < dim; ++k) {
        tempMat[k] = (double*)calloc(dim, sizeof(double));
        for (int l = 0; l < dim; ++l) {
            tempMat[k][l] = eigenVectors[l][k];
        }
    }
    for (int i = 0; i < dim; ++i) {
        for (int j = 0; j < dim; ++j) {
            eigenVectors[i][j] = tempMat[i][j];
        }
    }
    freeDMat(dim, tempMat);
    return eigenvalues;
}

_Complex double* cDiagJacobi(int dim, _Complex double** matrix, _Complex double** eigenVectors, double eps) {
    for (int i = 0; i < dim; ++i) {
        for (int j = 0; j < dim; ++j) {
            if (i == j) {
                eigenVectors[i][j] = 1. + 0.i;
            } else {
                eigenVectors[i][j] = 0. + 0.i;
            }
        }
    }
    double sn = 0.;
    for (int i = 0; i < dim; ++i) {
        for (int j = 0; j < dim; ++j) {
            if (i != j) {
                sn += cabs(matrix[i][j]);
            }
        }
    }
    sn /= (dim*dim);
    _Complex double** * cpyMat = cMatCpy(simDim(2, dim, dim), matrix);
    while (sn > eps) {
        int* loc = (int*)calloc(2, sizeof(int));
        loc[0] = 0, loc[1] = 1;
        for (int i = 0; i < dim-1; ++i) {
            for (int j = i+1; j < dim; ++j) {
                if (cabs((*cpyMat)[i][j]) > cabs((*cpyMat)[loc[0]][loc[1]])) {
                    loc[0] = i;
                    loc[1] = j;
                }
            }
        }
        double phase = .5 * carg((*cpyMat)[loc[0]][loc[1]]);
        double delta = pow((creal((*cpyMat)[loc[1]][loc[1]]) - creal((*cpyMat)[loc[0]][loc[0]]))/cabs((*cpyMat)[loc[0]][loc[1]]), 2.) + 4.;
        double tanT = (creal((*cpyMat)[loc[0]][loc[0]]) - creal((*cpyMat)[loc[1]][loc[1]]))/cabs((*cpyMat)[loc[0]][loc[1]]);
        if (fabs(tanT + sqrt(delta)) < fabs(tanT - sqrt(delta))) {
            tanT += sqrt(delta);
        } else {
            tanT -= sqrt(delta);
        }
        tanT /= 2.;
        double sinT = tanT / sqrt(1. + tanT*tanT);
        double cosT = 1. / sqrt(1. + tanT*tanT);
        _Complex double** * iterMat = cMatCpy(simDim(2, dim, dim), *cpyMat);
        for (int i = 0; i < dim; ++i) {
            (*cpyMat)[loc[0]][i] = (*iterMat)[loc[0]][i]*cosT*cexp((0. - 1.i)*phase) - (*iterMat)[loc[1]][i]*sinT*cexp((0. + 1.i)*phase);
            (*cpyMat)[loc[1]][i] = (*iterMat)[loc[1]][i]*cosT*cexp((0. + 1.i)*phase) + (*iterMat)[loc[0]][i]*sinT*cexp((0. - 1.i)*phase);
            (*cpyMat)[i][loc[0]] = (*iterMat)[i][loc[0]]*cosT*cexp((0. + 1.i)*phase) - (*iterMat)[i][loc[1]]*sinT*cexp((0. - 1.i)*phase);
            (*cpyMat)[i][loc[1]] = (*iterMat)[i][loc[1]]*cosT*cexp((0. - 1.i)*phase) + (*iterMat)[i][loc[0]]*sinT*cexp((0. + 1.i)*phase);
        }
        (*cpyMat)[loc[0]][loc[0]] = ((*iterMat)[loc[0]][loc[0]]*cosT*cexp((0. - 1.i)*phase) - (*iterMat)[loc[1]][loc[0]]*sinT*cexp((0. + 1.i)*phase))*cosT*cexp((0. + 1.i)*phase) -
                                 ((*iterMat)[loc[0]][loc[1]]*cosT*cexp((0. - 1.i)*phase) - (*iterMat)[loc[1]][loc[1]]*sinT*cexp((0. + 1.i)*phase))*sinT*cexp((0. - 1.i)*phase);
        (*cpyMat)[loc[0]][loc[1]] = ((*iterMat)[loc[0]][loc[1]]*cosT*cexp((0. - 1.i)*phase) - (*iterMat)[loc[1]][loc[1]]*sinT*cexp((0. + 1.i)*phase))*cosT*cexp((0. - 1.i)*phase) +
                                 ((*iterMat)[loc[0]][loc[0]]*cosT*cexp((0. - 1.i)*phase) - (*iterMat)[loc[1]][loc[0]]*sinT*cexp((0. + 1.i)*phase))*sinT*cexp((0. + 1.i)*phase);
        (*cpyMat)[loc[1]][loc[0]] = ((*iterMat)[loc[1]][loc[0]]*cosT*cexp((0. + 1.i)*phase) + (*iterMat)[loc[0]][loc[0]]*sinT*cexp((0. - 1.i)*phase))*cosT*cexp((0. + 1.i)*phase) -
                                 ((*iterMat)[loc[1]][loc[1]]*cosT*cexp((0. + 1.i)*phase) + (*iterMat)[loc[0]][loc[1]]*sinT*cexp((0. - 1.i)*phase))*sinT*cexp((0. - 1.i)*phase);
        (*cpyMat)[loc[1]][loc[1]] = ((*iterMat)[loc[1]][loc[1]]*cosT*cexp((0. + 1.i)*phase) + (*iterMat)[loc[0]][loc[1]]*sinT*cexp((0. - 1.i)*phase))*cosT*cexp((0. - 1.i)*phase) +
                                 ((*iterMat)[loc[1]][loc[0]]*cosT*cexp((0. + 1.i)*phase) + (*iterMat)[loc[0]][loc[0]]*sinT*cexp((0. - 1.i)*phase))*sinT*cexp((0. + 1.i)*phase);
        freeCMat(dim, *iterMat);
        sn = 0.;
        for (int i = 0; i < dim; ++i) {
            for (int j = 0; j < dim; ++j) {
                if (i != j) {
                    sn += cabs((*cpyMat)[i][j]);
                }
            }
        }
        sn /= (dim*dim);
        iterMat = cMatCpy(simDim(2, dim, dim), eigenVectors);
        for (int i = 0; i < dim; ++i) {
            eigenVectors[i][loc[0]] = (*iterMat)[i][loc[0]]*cosT*cexp((0. + 1.i)*phase) - (*iterMat)[i][loc[1]]*sinT*cexp((0. - 1.i)*phase);
            eigenVectors[i][loc[1]] = (*iterMat)[i][loc[1]]*cosT*cexp((0. - 1.i)*phase) + (*iterMat)[i][loc[0]]*sinT*cexp((0. + 1.i)*phase);
        }
        freeCMat(dim, *iterMat);
        free(loc);
    }
    _Complex double* eigenvalues = (_Complex double*)calloc(dim, sizeof(_Complex double));
    for (int i = 0; i < dim; ++i) {
        eigenvalues[i] = creal((*cpyMat)[i][i]) + 0.i;
    }
    freeCMat(dim, *cpyMat);
    return eigenvalues;
}

void dSvd(int* dims, double** matrix, double* sigma, double** uuL, double** vvR, double eps, double*(*dMatDiagMathod)()) {
    srand(time(NULL));
    double** * stdShapeMat; int ld, sd, rank = 0;
    if (dims[1] > dims[0]) {
        stdShapeMat = dMatCpy(dims, matrix);
        ld = dims[1];
        sd = dims[0];
    } else {
        *stdShapeMat = dCvtM().trans(simDim(2, dims[0], dims[1]), matrix);
        ld = dims[0];
        sd = dims[1];
        double** temp = uuL;
        uuL = vvR;
        vvR = temp;
    }
    double** * symmetryMat = dMulN().inner(simDim(4, sd, ld, sd, ld), *stdShapeMat, *stdShapeMat);
    double* lambdas = dMatDiagMathod(ld, *symmetryMat, vvR, eps);
    freeDMat(ld, *symmetryMat);
    for (int i = 0; i < ld-1; ++i) {
        for (int j = i+1; j < ld; ++j) {
            if (lambdas[i] < lambdas[j]) {
                double temp = lambdas[i];
                lambdas[i] = lambdas[j];
                lambdas[j] = temp;
                dSwtM('R').exc(vvR, simDim(2, ld, ld), simDim(2, i, j));
            }
        }
    }
    for (int i = 0; i < sd; ++i) {
        if (lambdas[i] > eps) {
            sigma[i] = sqrt(lambdas[i]);
            rank++;
        } else {
            sigma[i] = 0.;
        }
    }
    freeDVec(lambdas);
    double** * tu = dMulN().multiple(simDim(4, sd, ld, ld, ld), *stdShapeMat, vvR);
    for (int i = 0; i < sd; ++i) {
        for (int j = 0; j < sd; ++j) {
            uuL[j][i] = *tu[j][i] / sigma[i];
        }
    }
    if (rank < sd) {
        for (int i = rank; i < sd; ++i) {
            for (int j = 0; j < sd; ++j) {
                uuL[j][i] = 1.*rand()/RAND_MAX;
            }
        }
        dSchmidt(sd, uuL);
    }
    if (dims[1] < dims[0]) {
        double** temp = uuL;
        uuL = vvR;
        vvR = temp;
    }
    freeDMat(sd, *stdShapeMat);
    freeDMat(sd, *tu);
}

void cSvd(int* dims, _Complex double** matrix, double* sigma, _Complex double** uuL, _Complex double** vvR, double eps, _Complex double*(*cMatDiagMathod)()) {
    srand(time(NULL));
    _Complex double** * stdShapeMat; int ld, sd, rank = 0;
    if (dims[1] > dims[0]) {
        stdShapeMat = cMatCpy(dims, matrix);
        ld = dims[1];
        sd = dims[0];
    } else {
        *stdShapeMat = cCvtM().dag(simDim(2, dims[0], dims[1]), matrix);
        ld = dims[0];
        sd = dims[1];
        _Complex double** temp = uuL;
        uuL = vvR;
        vvR = temp;
    }
    _Complex double** * hermiteMat = cMulN().inner(simDim(4, sd, ld, sd, ld), *stdShapeMat, *stdShapeMat);
    _Complex double* lambdas = cMatDiagMathod(ld, *hermiteMat, vvR, eps);
    freeCMat(ld, *hermiteMat);
    for (int i = 0; i < ld-1; ++i) {
        for (int j = i+1; j < ld; ++j) {
            if (cabs(lambdas[i]) < cabs(lambdas[j])) {
                _Complex double temp = lambdas[i];
                lambdas[i] = lambdas[j];
                lambdas[j] = temp;
                cSwtM('R').exc(vvR, simDim(2, ld, ld), simDim(2, i, j));
            }
        }
    }
    for (int i = 0; i < sd; ++i) {
        if (cabs(lambdas[i]) > eps) {
            sigma[i] = sqrt(lambdas[i]);
            rank++;
        } else {
            sigma[i] = 0.;
        }
    }
    freeCVec(lambdas);
    _Complex double** * tu = cMulN().multiple(simDim(4, sd, ld, ld, ld), *stdShapeMat, vvR);
    for (int i = 0; i < sd; ++i) {
        for (int j = 0; j < sd; ++j) {
            uuL[j][i] = *tu[j][i]/(sigma[i] + 0.i);
        }
    }
    if (rank < sd) {
        for (int i = rank; i < sd; ++i) {
            for (int j = 0; j < sd; ++j) {
                uuL[j][i] = 1.*rand()/RAND_MAX + 1.i*rand()/RAND_MAX;
            }
        }
        cSchmidt(sd, uuL);
    }
    if (dims[1] < dims[0]) {
        _Complex double** temp = uuL;
        uuL = vvR;
        vvR = temp;
    }
    freeCMat(sd, *stdShapeMat);
    freeCMat(sd, *tu);
}

double* dDiagPower(int dim, double** matrix, double** eigenVectors, double eps) {
    srand(time(NULL));
    double* eigenvalues = (double*)calloc(1, sizeof(double));
    double temp = 0.;
    for (int i = 0; i < dim; ++i) {
        eigenVectors[i][0] = 1.*rand();
        temp += eigenVectors[i][0]*eigenVectors[i][0];
    }
    for (int i = 0; i < dim; ++i) {
        eigenVectors[i][0] /= sqrt(temp);
    }
    double* tempVec = (double*)calloc(dim, sizeof(double));
    double delta = 1.e+24;
    do {
        temp = 0.;
        for (int i = 0; i < dim; ++i) {
            tempVec[i] = 0.;
            for (int j = 0; j < dim; ++j) {
                tempVec[i] += matrix[i][j]*eigenVectors[j][0];
            }
            temp += tempVec[i]*tempVec[i];
        }
        delta = eigenvalues[0];
        eigenvalues[0] = 0.;
        for (int i = 0; i < dim; ++i) {
            eigenVectors[i][0] = tempVec[i]/sqrt(temp);
            for (int j = 0; j < dim; ++j) {
                eigenvalues[0] += tempVec[i]/sqrt(temp)*matrix[i][j]*tempVec[j]/sqrt(temp);
            }
        }
        delta -= eigenvalues[0];
    } while (fabs(delta) > eps);
    freeDVec(tempVec);
    return eigenvalues;
}

_Complex double* cDiagPower(int dim, _Complex double** matrix, _Complex double** eigenVectors, double eps) {
    srand(time(NULL));
    _Complex double* eigenvalues = (_Complex double*)calloc(1, sizeof(_Complex double));
    for (int i = 0; i < dim; ++i) {
        eigenVectors[i][0] = 1.*rand() + 1.i*rand();
    }
    _Complex double* tempVec = (_Complex double*)calloc(dim, sizeof(_Complex double));
    _Complex double temp = 0. + 0.i;
    _Complex double delta = 1.e+24 - 1.e-9i;
    do {
        temp = 0. + 0.i;
        for (int i = 0; i < dim; ++i) {
            tempVec[i] = 0.;
            for (int j = 0; j < dim; ++j) {
                tempVec[i] += matrix[i][j]*eigenVectors[j][0];
            }
            temp += conj(tempVec[i])*tempVec[i];
        }
        delta = eigenvalues[0];
        eigenvalues[0] = 0. + 0.i;
        for (int i = 0; i < dim; ++i) {
            eigenVectors[i][0] = tempVec[i]/sqrt(creal(temp));
            for (int j = 0; j < dim; ++j) {
                eigenvalues[0] += conj(tempVec[i]/sqrt(creal(temp)))*matrix[i][j]*tempVec[j]/sqrt(creal(temp));
            }
        }
        delta -= eigenvalues[0];
    } while (cabs(delta) > eps);
    freeCVec(tempVec);
    return eigenvalues;
}

double* dDiagPowerWithShift(int dim, double** matrix, double** eigenVectors, double shift, double eps) {
    srand(time(NULL));
    double* eigenvalues = (double*)calloc(1, sizeof(double));
    double temp = 0.;
    for (int i = 0; i < dim; ++i) {
        eigenVectors[i][0] = 1.*rand();
        temp += eigenVectors[i][0]*eigenVectors[i][0];
    }
    for (int i = 0; i < dim; ++i) {
        eigenVectors[i][0] /= sqrt(temp);
    }
    double* tempVec = (double*)calloc(dim, sizeof(double));
    double delta = 1.e+24;
    do {
        temp = 0.;
        for (int i = 0; i < dim; ++i) {
            tempVec[i] = 0.;
            for (int j = 0; j < dim; ++j) {
                tempVec[i] += (matrix[i][j] - shift*(i == j))*eigenVectors[j][0];
            }
            temp += tempVec[i]*tempVec[i];
        }
        delta = eigenvalues[0];
        eigenvalues[0] = 0.;
        for (int i = 0; i < dim; ++i) {
            eigenVectors[i][0] = tempVec[i]/sqrt(temp);
            for (int j = 0; j < dim; ++j) {
                eigenvalues[0] += tempVec[i]/sqrt(temp)*(matrix[i][j] - shift*(i == j))*tempVec[j]/sqrt(temp);
            }
        }
        delta -= eigenvalues[0];
    } while (fabs(delta) > eps);
    freeDVec(tempVec);
    eigenvalues[0] += shift;
    return eigenvalues;
}

_Complex double* cDiagPowerWithShift(int dim, _Complex double** matrix, _Complex double** eigenVectors, _Complex double shift, double eps) {
    srand(time(NULL));
    _Complex double* eigenvalues = (_Complex double*)calloc(1, sizeof(_Complex double));
    for (int i = 0; i < dim; ++i) {
        eigenVectors[i][0] = 1.*rand() + 1.i*rand();
    }
    _Complex double* tempVec = (_Complex double*)calloc(dim, sizeof(_Complex double));
    _Complex double temp = 0. + 0.i;
    _Complex double delta = 1.e+24 - 1.e+9i;
    do {
        temp = 0. + 0.i;
        for (int i = 0; i < dim; ++i) {
            tempVec[i] = 0.;
            for (int j = 0; j < dim; ++j) {
                tempVec[i] += (matrix[i][j] - shift*(i == j))*eigenVectors[j][0];
            }
            temp += conj(tempVec[i])*tempVec[i];
        }
        delta = eigenvalues[0];
        eigenvalues[0] = 0. + 0.i;
        for (int i = 0; i < dim; ++i) {
            eigenVectors[i][0] = tempVec[i]/sqrt(creal(temp));
            for (int j = 0; j < dim; ++j) {
                eigenvalues[0] += conj(tempVec[i]/sqrt(creal(temp)))*(matrix[i][j] - shift*(i == j))*tempVec[j]/sqrt(creal(temp));
            }
        }
        delta -= eigenvalues[0];
    } while (cabs(delta) > eps);
    freeCVec(tempVec);
    eigenvalues[0] += shift;
    return eigenvalues;
}

double* dDiagAntiPowerWithShift(int dim, double** matrix, double** eigenVectors, double scalar, double eps) {
    srand(time(NULL));
    double* eigenvalues = (double*)calloc(1, sizeof(double));
    double temp = 0.;
    for (int i = 0; i < dim; ++i) {
        eigenVectors[i][0] = 1.*rand();
        temp += eigenVectors[i][0]*eigenVectors[i][0];
    }
    for (int i = 0; i < dim; ++i) {
        eigenVectors[i][0] /= sqrt(temp);
    }
    double** tempMat = dNewM(SHAPE_RECTANGLE, 0).create_rm(simDim(2, dim, dim+1));
    for (int i = 0; i < dim; ++i) {
        for (int j = 0; j < dim; ++j) {
            tempMat[i][j] = matrix[i][j] - scalar*(i == j);
        }
    }
    double* tempVec = (double*)calloc(dim, sizeof(double));
    double delta = 1.e+24;
    do {
        temp = 0.;
        for (int i = 0; i < dim; ++i) {
            tempMat[i][dim] = eigenVectors[i][0];
        }
        double** solve = dMatrixLineSimBeta(simDim(2, dim, dim+1), tempMat);
        for (int i = 0; i < dim; ++i) {
            tempVec[i] = solve[i][dim];
            temp += tempVec[i]*tempVec[i];
        }
        freeDMat(dim, solve);
        delta = eigenvalues[0];
        eigenvalues[0] = 0.;
        for (int i = 0; i < dim; ++i) {
            eigenVectors[i][0] = tempVec[i]/sqrt(temp);
            for (int j = 0; j < dim; ++j) {
                eigenvalues[0] += tempVec[i]/sqrt(temp)*matrix[i][j]*tempVec[j]/sqrt(temp);
            }
        }
        delta -= eigenvalues[0];
    } while (fabs(delta) > eps);
    freeDVec(tempVec);
    freeDMat(dim, tempMat);
    return eigenvalues;
}

_Complex double* cDiagAntiPowerWithShift(int dim, _Complex double** matrix, _Complex double** eigenVectors, _Complex double scalar, double eps) {
    srand(time(NULL));
    _Complex double* eigenvalues = (_Complex double*)calloc(1, sizeof(_Complex double));
    _Complex double temp = 0. + 0.i;
    for (int i = 0; i < dim; ++i) {
        eigenVectors[i][0] = 1.*rand() + 1.i*rand();
        temp += conj(eigenVectors[i][0])*eigenVectors[i][0];
    }
    for (int i = 0; i < dim; ++i) {
        eigenVectors[i][0] /= sqrt(creal(temp));
    }
    _Complex double** tempMat = cNewM(SHAPE_RECTANGLE, 0).create_rm(simDim(2, dim, dim+1));
    for (int i = 0; i < dim; ++i) {
        for (int j = 0; j < dim; ++j) {
            tempMat[i][j] = matrix[i][j] - scalar*(i == j);
        }
    }
    _Complex double* tempVec = (_Complex double*)calloc(dim, sizeof(_Complex double));
    _Complex double delta = 1.e+24 - 1.e+9i;
    do {
        temp = 0.;
        for (int i = 0; i < dim; ++i) {
            tempMat[i][dim] = eigenVectors[i][0];
        }
        _Complex double** solve = cMatrixLineSimBeta(simDim(2, dim, dim+1), tempMat);
        for (int i = 0; i < dim; ++i) {
            tempVec[i] = solve[i][dim];
            temp += conj(tempVec[i])*tempVec[i];
        }
        freeCMat(dim, solve);
        delta = eigenvalues[0];
        eigenvalues[0] = 0.;
        for (int i = 0; i < dim; ++i) {
            eigenVectors[i][0] = tempVec[i]/sqrt(creal(temp));
            for (int j = 0; j < dim; ++j) {
                eigenvalues[0] += conj(tempVec[i]/sqrt(creal(temp)))*matrix[i][j]*tempVec[j]/sqrt(creal(temp));
            }
        }
        delta -= eigenvalues[0];
    } while (cabs(delta) > eps);
    freeCVec(tempVec);
    freeCMat(dim, tempMat);
    return eigenvalues;
}

double* dDiagRayleighQI(int dim, double** matrix, double** eigenVectors, double eps) {
    srand(time(NULL));
    double* eigenvalues = (double*)calloc(1, sizeof(double));
    double temp = 0.;
    for (int i = 0; i < dim; ++i) {
        eigenVectors[i][0] = 1.*rand();
        temp += eigenVectors[i][0]*eigenVectors[i][0];
    }
    for (int i = 0; i < dim; ++i) {
        eigenVectors[i][0] /= sqrt(temp);
    }
    double sigma = 0.;
    for (int i = 0; i < dim; ++i) {
        for (int j = 0; j < dim; ++j) {
            sigma += eigenVectors[i][0]*matrix[i][j]*eigenVectors[j][0];
        }
    }
    double** tempMat = dNewM(SHAPE_RECTANGLE, 0).create_rm(simDim(2, dim, dim+1));
    for (int i = 0; i < dim; ++i) {
        for (int j = 0; j < dim; ++j) {
            tempMat[i][j] = matrix[i][j] - sigma*(i == j);
        }
    }
    double* tempVec = (double*)calloc(dim, sizeof(double));
    eigenvalues[0] = sigma;
    do {
        temp = 0.;
        for (int i = 0; i < dim; ++i) {
            tempMat[i][i] += (eigenvalues[0] - sigma);
            tempMat[i][dim] = eigenVectors[i][0];
        }
        double** solve = dMatrixLineSimBeta(simDim(2, dim, dim+1), tempMat);
        for (int i = 0; i < dim; ++i) {
            tempVec[i] = solve[i][dim];
            temp += tempVec[i]*tempVec[i];
        }
        freeDMat(dim, solve);
        eigenvalues[0] = sigma;
        sigma = 0.;
        for (int i = 0; i < dim; ++i) {
            eigenVectors[i][0] = tempVec[i]/sqrt(temp);
            for (int j = 0; j < dim; ++j) {
                sigma += tempVec[i]/sqrt(temp)*matrix[i][j]*tempVec[j]/sqrt(temp);
            }
        }
    } while (fabs(sigma - eigenvalues[0]) > eps);
    freeDVec(tempVec);
    freeDMat(dim, tempMat);
    return eigenvalues;
}

_Complex double* cDiagRayleighQI(int dim, _Complex double** matrix, _Complex double** eigenVectors, double eps) {
    srand(time(NULL));
    _Complex double* eigenvalues = (_Complex double*)calloc(1, sizeof(_Complex double));
    _Complex double temp = 0. + 0.i;
    for (int i = 0; i < dim; ++i) {
        eigenVectors[i][0] = 1.*rand() + 1.i*rand();
        temp += conj(eigenVectors[i][0])*eigenVectors[i][0];
    }
    for (int i = 0; i < dim; ++i) {
        eigenVectors[i][0] /= sqrt(creal(temp));
    }
    _Complex double sigma = 0. + 0.i;
    for (int i = 0; i < dim; ++i) {
        for (int j = 0; j < dim; ++j) {
            sigma += conj(eigenVectors[i][0])*matrix[i][j]*eigenVectors[j][0];
        }
    }
    _Complex double** tempMat = cNewM(SHAPE_RECTANGLE, 0).create_rm(simDim(2, dim, dim+1));
    for (int i = 0; i < dim; ++i) {
        for (int j = 0; j < dim; ++j) {
            tempMat[i][j] = matrix[i][j] - sigma*(i == j);
        }
    }
    _Complex double* tempVec = (_Complex double*)calloc(dim, sizeof(_Complex double));
    eigenvalues[0] = sigma;
    do {
        temp = 0. + 0.i;
        for (int i = 0; i < dim; ++i) {
            tempMat[i][i] += (eigenvalues[0] - sigma);
            tempMat[i][dim] = eigenVectors[i][0];
        }
        _Complex double** solve = cMatrixLineSimBeta(simDim(2, dim, dim+1), tempMat);
        for (int i = 0; i < dim; ++i) {
            tempVec[i] = solve[i][dim];
            temp += conj(tempVec[i])*tempVec[i];
        }
        freeCMat(dim, solve);
        eigenvalues[0] = sigma;
        sigma = 0. + 0.i;
        for (int i = 0; i < dim; ++i) {
            eigenVectors[i][0] = tempVec[i]/sqrt(creal(temp));
            for (int j = 0; j < dim; ++j) {
                sigma += conj(tempVec[i])/sqrt(creal(temp))*matrix[i][j]*tempVec[j]/sqrt(creal(temp));
            }
        }
    } while (cabs(sigma - eigenvalues[0]) > eps);
    freeCVec(tempVec);
    freeCMat(dim, tempMat);
    return eigenvalues;
}

double* dDiagQR(int dim, double** matrix, double** eigenVectors, double eps) {
    double* eigenvalues = (double*)calloc(dim, sizeof(double));
    for (int i = 0; i < dim; ++i) {
        for (int j = 0; j < dim; ++j) {
            eigenVectors[i][j] = 1.*(i == j);
        }
    }
    double** * matrixK;
    matrixK = dMatCpy(simDim(2, dim, dim), matrix);
    double delta = 1.e+24;
    do {
        double** * matrixQ = dMatCpy(simDim(2, dim, dim), *matrixK);
        dSchmidt(dim, *matrixQ);
        double** * tempMat = dMulN().multiple(simDim(4, dim, dim, dim, dim), *dMulN().inner(simDim(4, dim, dim, dim, dim), *matrixQ, *matrixK), *matrixQ);
        delta = 0.;
        for (int i = 0; i < dim; ++i) {
            for (int j = 0; j < dim; ++j) {
                (*matrixK)[i][j] = (*tempMat)[i][j];
                (*tempMat)[i][j] = eigenVectors[i][j];
            }
            delta += fabs(eigenvalues[i] - (*matrixK)[i][i]);
            eigenvalues[i] = (*matrixK)[i][i];
        }
        for (int i = 0; i < dim; ++i) {
            for (int j = 0; j < dim; ++j) {
                eigenVectors[i][j] = 0.;
                for (int k = 0; k < dim; ++k) {
                    eigenVectors[i][j] += (*tempMat)[i][k]*(*matrixQ)[k][j];
                }
            }
        }
        freeDMat(dim, *matrixQ);
        freeDMat(dim, *tempMat);
    } while (delta/dim > eps);
    freeDMat(dim, *matrixK);
    return eigenvalues;
}

_Complex double* cDiagQR(int dim, _Complex double** matrix, _Complex double** eigenVectors, double eps) {
    _Complex double* eigenvalues = (_Complex double*)calloc(dim, sizeof(_Complex double));
    for (int i = 0; i < dim; ++i) {
        for (int j = 0; j < dim; ++j) {
            eigenVectors[i][j] = 1.*(i == j) + 0.i;
        }
    }
    _Complex double** * matrixK;
    matrixK = cMatCpy(simDim(2, dim, dim), matrix);
    _Complex double delta = 1.e+24 - 1.e+9i;
    do {
        _Complex double** * matrixQ = cMatCpy(simDim(2, dim, dim), *matrixK);
        cSchmidt(dim, *matrixQ);
        _Complex double** * tempMat = cMulN().multiple(simDim(4, dim, dim, dim, dim), *cMulN().inner(simDim(4, dim, dim, dim, dim), *matrixQ, *matrixK), *matrixQ);
        delta = 0. + 0.i;
        for (int i = 0; i < dim; ++i) {
            for (int j = 0; j < dim; ++j) {
                (*matrixK)[i][j] = (*tempMat)[i][j];
                (*tempMat)[i][j] = eigenVectors[i][j];
            }
            delta += cabs(eigenvalues[i] - (*matrixK)[i][i]);
            eigenvalues[i] = (*matrixK)[i][i];
        }
        for (int i = 0; i < dim; ++i) {
            for (int j = 0; j < dim; ++j) {
                eigenVectors[i][j] = 0. + 0.i;
                for (int k = 0; k < dim; ++k) {
                    eigenVectors[i][j] += (*tempMat)[i][k]*(*matrixQ)[k][j];
                }
            }
        }
        freeCMat(dim, *matrixQ);
        freeCMat(dim, *tempMat);
    } while (cabs(delta)/dim > eps);
    freeCMat(dim, *matrixK);
    return eigenvalues;
}

double* dDiagQRWithAccAlpha(int dim, double** matrix, double** eigenVectors, double eps) {
    double* eigenvalues = (double*)calloc(dim, sizeof(double));
    for (int i = 0; i < dim; ++i) {
        for (int j = 0; j < dim; ++j) {
            eigenVectors[i][j] = 1.*(i == j);
        }
    }
    double** * matrixK;
    matrixK = dMatCpy(simDim(2, dim, dim), matrix);
    double delta = 1.e+24;
    do {
        double** * matrixQ = dMatCpy(simDim(2, dim, dim), *matrixK);
        double shift = (*matrixQ)[dim-1][dim-1];
        for (int i = 0; i < dim; ++i) {
            (*matrixQ)[i][i] -= shift;
            (*matrixK)[i][i] -= shift;
        }
        dSchmidt(dim, *matrixQ);
        double** * tempMat = dMulN().multiple(simDim(4, dim, dim, dim, dim), *dMulN().inner(simDim(4, dim, dim, dim, dim), *matrixQ, *matrixK), *matrixQ);
        for (int i = 0; i < dim; ++i) {
            (*tempMat)[i][i] += shift;
        }
        delta = 0.;
        for (int i = 0; i < dim; ++i) {
            for (int j = 0; j < dim; ++j) {
                (*matrixK)[i][j] = (*tempMat)[i][j];
                (*tempMat)[i][j] = eigenVectors[i][j];
            }
            delta += fabs(eigenvalues[i] - (*matrixK)[i][i]);
            eigenvalues[i] = (*matrixK)[i][i];
        }
        for (int i = 0; i < dim; ++i) {
            for (int j = 0; j < dim; ++j) {
                eigenVectors[i][j] = 0.;
                for (int k = 0; k < dim; ++k) {
                    eigenVectors[i][j] += (*tempMat)[i][k]*(*matrixQ)[k][j];
                }
            }
        }
        freeDMat(dim, *matrixQ);
        freeDMat(dim, *tempMat);
    } while (delta/dim > eps);
    freeDMat(dim, *matrixK);
    return eigenvalues;
}

_Complex double* cDiagQRWithAccAlpha(int dim, _Complex double** matrix, _Complex double** eigenVectors, double eps) {
    _Complex double* eigenvalues = (_Complex double*)calloc(dim, sizeof(_Complex double));
    for (int i = 0; i < dim; ++i) {
        for (int j = 0; j < dim; ++j) {
            eigenVectors[i][j] = 1.*(i == j);
        }
    }
    _Complex double** * matrixK;
    matrixK = cMatCpy(simDim(2, dim, dim), matrix);
    _Complex double delta = 1.e+24 - 1.e+9i;
    do {
        _Complex double** * matrixQ = cMatCpy(simDim(2, dim, dim), *matrixK);
        _Complex double shift = (*matrixQ)[dim-1][dim-1];
        for (int i = 0; i < dim; ++i) {
            (*matrixQ)[i][i] -= shift;
            (*matrixK)[i][i] -= shift;
        }
        cSchmidt(dim, *matrixQ);
        _Complex double** * tempMat = cMulN().multiple(simDim(4, dim, dim, dim, dim), *cMulN().inner(simDim(4, dim, dim, dim, dim), *matrixQ, *matrixK), *matrixQ);
        for (int i = 0; i < dim; ++i) {
            (*tempMat)[i][i] += shift;
        }
        delta = 0. + 0.i;
        for (int i = 0; i < dim; ++i) {
            for (int j = 0; j < dim; ++j) {
                (*matrixK)[i][j] = (*tempMat)[i][j];
                (*tempMat)[i][j] = eigenVectors[i][j];
            }
            delta += cabs(eigenvalues[i] - (*matrixK)[i][i]);
            eigenvalues[i] = (*matrixK)[i][i];
        }
        for (int i = 0; i < dim; ++i) {
            for (int j = 0; j < dim; ++j) {
                eigenVectors[i][j] = 0.;
                for (int k = 0; k < dim; ++k) {
                    eigenVectors[i][j] += (*tempMat)[i][k]*(*matrixQ)[k][j];
                }
            }
        }
        freeCMat(dim, *matrixQ);
        freeCMat(dim, *tempMat);
    } while (cabs(delta)/dim > eps);
    freeCMat(dim, *matrixK);
    return eigenvalues;
}

double* dDiagQRWithAccBeta(int dim, double** matrix, double** eigenVectors, double eps) {
    srand(time(NULL));
    double** objVal = *dMatCpy(simDim(2, dim, dim), matrix);
    double** objVec = dNewM(SHAPE_SQUARE, 'U').create_sm(dim);
    double** tempMat = dNewM(SHAPE_RECTANGLE, 0).create_rm(simDim(2, dim, dim));
    for (int i = 0; i < dim-2; ++i) {
        double* tmpHouseV = (double*)calloc(dim-i-1, sizeof(double));
        for (int j = 0; j < dim-i-1; ++j) {
            tmpHouseV[j] = objVal[i+j+1][i];
        }
        for (int j = 0; j < dim-i-1; ++j) {
            for (int k = 0; k < dim-i; ++k) {
                tempMat[j][k] = objVal[j+i+1][k+i];
            }
        }
        for (int j = 0; j < dim-i-1; ++j) {
            for (int k = 0; k < dim-i; ++k) {
                double tempSum = 0.;
                for (int l = 0; l < dim-i-1; ++l) {
                    tempSum += tmpHouseV[l]*tempMat[l][k];
                }
                objVal[j+i+1][k+i] -= 2.*tmpHouseV[j]*tempSum/dLenVec(dim-i-1, tmpHouseV)/dLenVec(dim-i-1, tmpHouseV);
            }
        }
        for (int j = 0; j < dim; ++j) {
            for (int k = 0; k < dim-i-1; ++k) {
                tempMat[j][k] = objVal[j][k+i+1];
            }
        }
        for (int j = 0; j < dim; ++j) {
            for (int k = 0; k < dim-i-1; ++k) {
                double tempSum = 0.;
                for (int l = 0; l < dim-i-1; ++l) {
                    tempSum += tempMat[j][l]*tmpHouseV[l];
                }
                objVal[j][k+i+1] -= 2.*tempSum*tmpHouseV[k]/dLenVec(dim-i-1, tmpHouseV)/dLenVec(dim-i-1, tmpHouseV);
            }
        }
        for (int j = 0; j < dim-i-1; ++j) {
            for (int k = 0; k < dim-i; ++k) {
                tempMat[j][k] = objVec[j+i+1][k+i];
            }
        }
        for (int j = 0; j < dim-i-1; ++j) {
            for (int k = 0; k < dim-i; ++k) {
                double tempSum = 0.;
                for (int l = 0; l < dim-i-1; ++l) {
                    tempSum += tmpHouseV[l]*tempMat[l][k];
                }
                objVec[j+i+1][k+i] -= 2.*tmpHouseV[j]*tempSum/dLenVec(dim-i-1, tmpHouseV)/dLenVec(dim-i-1, tmpHouseV);
            }
        }
        freeDVec(tmpHouseV);
    }
    double delta = 1.e+24;
    double* retVals = dNewV(0).create(dim);
    do {
        for (int i = 0; i < dim; ++i) {
            retVals[i] = objVal[i][i];
        }
        double shift = objVal[dim-1][dim-1];
        double** mh = dNewM(SHAPE_RECTANGLE, 0).create_rm(simDim(2, 3, 3));
        mh[0][0] = objVal[0][0]*objVal[0][0] + objVal[0][1]*objVal[1][0] - 2.*shift*objVal[0][0] + pow(fabs(shift), 2.);
        mh[1][0] = objVal[1][0]*(objVal[0][0] + objVal[1][1] - 2.*shift);
        mh[2][0] = objVal[1][0]*objVal[2][1];
        for (int i = 0; i < 3; ++i) {
            mh[i][0] = 1.*rand()/RAND_MAX;
            mh[i][1] = 1.*rand()/RAND_MAX;
        }
        dSchmidt(3, mh);
        for (int i = 0; i < dim; ++i) {
            tempMat[0][i] = objVal[0][i];
            tempMat[1][i] = objVal[1][i];
            tempMat[2][i] = objVal[2][i];
        }
        for (int i = 0; i < 3; ++i) {
            for (int j = 0; j < 3; ++j) {
                objVal[i][j] = 0.;
                for (int k = 0; k < 3; ++k) {
                    objVal[i][j] += mh[i][k]*tempMat[k][j];
                }
            }
            for (int j = 3; j < dim; ++j) {
                objVal[i][j] = 0.;
                for (int k = 0; k < 3; ++k) {
                    objVal[i][j] += mh[i][k]*tempMat[k][j];
                }
            }
        }
        for (int i = 0; i < dim; ++i) {
            tempMat[i][0] = objVal[i][0];
            tempMat[i][1] = objVal[i][1];
            tempMat[i][2] = objVal[i][2];
        }
        for (int i = 0; i < 3; ++i) {
            for (int j = 0; j < 3; ++j) {
                objVal[i][j] = 0.;
                for (int k = 0; k < 3; ++k) {
                    objVal[i][j] += tempMat[i][k]*mh[j][k];
                }
            }
        }
        for (int i = 3; i < dim; ++i) {
            for (int j = 0; j < 3; ++j) {
                objVal[i][j] = 0.;
                for (int k = 0; k < 3; ++k) {
                    objVal[i][j] += tempMat[i][k]*mh[j][k];
                }
            }
        }
        for (int i = 0; i < dim; ++i) {
            tempMat[i][0] = objVec[i][0];
            tempMat[i][1] = objVec[i][1];
            tempMat[i][2] = objVec[i][2];
        }
        for (int i = 0; i < 3; ++i) {
            for (int j = 0; j < 3; ++j) {
                objVec[i][j] = 0.;
                for (int k = 0; k < 3; ++k) {
                    objVec[i][j] += tempMat[i][k]*mh[j][k];
                }
            }
        }
        for (int i = 3; i < dim; ++i) {
            for (int j = 0; j < 3; ++j) {
                objVec[i][j] = 0.;
                for (int k = 0; k < 3; ++k) {
                    objVec[i][j] += tempMat[i][k]*mh[j][k];
                }
            }
        }
        for (int i = 0; i < dim-3; ++i) {
            double* tmpHouseV = simDVec(3, objVal[i+1][i], objVal[i+2][i], objVal[i+3][i]);
            for (int j = i+1; j < i+4; ++j) {
                for (int k = 0; k < dim; ++k) {
                    tempMat[j][k] = objVal[j][k];
                }
            }
            for (int j = i+1; j < i+4; ++j) {
                for (int k = 0; k < i+1; ++k) {
                    double tempSum = 0.;
                    for (int l = 0; l < 3; ++l) {
                        tempSum += tmpHouseV[l]*tempMat[l+i+1][k];
                    }
                    objVal[j][k] -= 2.*tmpHouseV[j-i-1]*tempSum/dLenVec(3, tmpHouseV)/dLenVec(3, tmpHouseV);
                }
                for (int k = i+1; k < i+4; ++k) {
                    double tempSum = 0.;
                    for (int l = 0; l < 3; ++l) {
                        tempSum += tmpHouseV[l]*tempMat[l+i+1][k];
                    }
                    objVal[j][k] -= 2.*tmpHouseV[j-i-1]*tempSum/dLenVec(3, tmpHouseV)/dLenVec(3, tmpHouseV);
                }
                for (int k = i+4; k < dim; ++k) {
                    double tempSum = 0.;
                    for (int l = 0; l < 3; ++l) {
                        tempSum += tmpHouseV[l]*tempMat[l+i+1][k];
                    }
                    objVal[j][k] -= 2.*tmpHouseV[j-i-1]*tempSum/dLenVec(3, tmpHouseV)/dLenVec(3, tmpHouseV);
                }
            }
            for (int j = 0; j < dim; ++j) {
                for (int k = i+1; k < i+4; ++k) {
                    tempMat[j][k] = objVal[j][k];
                }
            }
            for (int j = 0; j < i+1; ++j) {
                for (int k = i+1; k < i+4; ++k) {
                    double tempSum = 0.;
                    for (int l = 0; l < 3; ++l) {
                        tempSum += objVal[j][l+i+1]*tmpHouseV[l];
                    }
                    objVal[j][k] -= 2.*tempSum*tmpHouseV[j-i-1]/dLenVec(3, tmpHouseV)/dLenVec(3, tmpHouseV);
                }
            }
            for (int j = i+1; j < i+4; ++j) {
                for (int k = i+1; k < i+4; ++k) {
                    double tempSum = 0.;
                    for (int l = 0; l < 3; ++l) {
                        tempSum += objVal[j][l+i+1]*tmpHouseV[l];
                    }
                    objVal[j][k] -= 2.*tempSum*tmpHouseV[j-i-1]/dLenVec(3, tmpHouseV)/dLenVec(3, tmpHouseV);
                }
            }
            for (int j = i+4; j < dim; ++j) {
                for (int k = i+1; k < i+4; ++k) {
                    double tempSum = 0.;
                    for (int l = 0; l < 3; ++l) {
                        tempSum += objVal[j][l+i+1]*tmpHouseV[l];
                    }
                    objVal[j][k] -= 2.*tempSum*tmpHouseV[j-i-1]/dLenVec(3, tmpHouseV)/dLenVec(3, tmpHouseV);
                }
            }
            for (int j = 0; j < dim; ++j) {
                for (int k = i+1; k < i+4; ++k) {
                    tempMat[j][k] = objVec[j][k];
                }
            }
            for (int j = 0; j < i+1; ++j) {
                for (int k = i+1; k < i+4; ++k) {
                    double tempSum = 0.;
                    for (int l = 0; l < 3; ++l) {
                        tempSum += objVec[j][l+i+1]*tmpHouseV[l];
                    }
                    objVec[j][k] -= 2.*tempSum*tmpHouseV[j-i-1]/dLenVec(3, tmpHouseV)/dLenVec(3, tmpHouseV);
                }
            }
            for (int j = i+1; j < i+4; ++j) {
                for (int k = i+1; k < i+4; ++k) {
                    double tempSum = 0.;
                    for (int l = 0; l < 3; ++l) {
                        tempSum += objVec[j][l+i+1]*tmpHouseV[l];
                    }
                    objVec[j][k] -= 2.*tempSum*tmpHouseV[j-i-1]/dLenVec(3, tmpHouseV)/dLenVec(3, tmpHouseV);
                }
            }
            for (int j = i+4; j < dim; ++j) {
                for (int k = i+1; k < i+4; ++k) {
                    double tempSum = 0.;
                    for (int l = 0; l < 3; ++l) {
                        tempSum += objVec[j][l+i+1]*tmpHouseV[l];
                    }
                    objVec[j][k] -= 2.*tempSum*tmpHouseV[j-i-1]/dLenVec(3, tmpHouseV)/dLenVec(3, tmpHouseV);
                }
            }
            freeDVec(tmpHouseV);
        }
        double*** giv = dCreateGivensMatrix(dim, simDim(2, dim-2, dim-1), objVal[dim-2][dim-2]);
        dGivensTransDirect(dim, &objVal, dCvtM().trans(simDim(2, dim, dim), *giv));
        for (int i = 0; i < dim; ++i) {
            for (int j = 0; j < dim; ++j) {
                tempMat[i][j] = objVec[i][j];
            }
        }
        for (int i = 0; i < dim; ++i) {
            for (int j = 0; j < dim; ++j) {
                objVec[i][j] = 0.;
                for (int k = 0; k < dim; ++k) {
                    objVec[i][j] += tempMat[i][k]*(*giv)[j][k];
                }
            }
        }
        freeDMat(dim, *giv);
        freeDMat(3, mh);
        delta = 0.;
        for (int i = 0; i < dim; ++i) {
            delta += fabs(retVals[i] - objVal[i][i]);
        }
    } while (delta/dim > eps);
    freeDMat(dim, tempMat);
    return retVals;
}

_Complex double* cDiagQRWithAccBeta(int dim, _Complex double** matrix, _Complex double** eigenVectors, double eps) {
    srand(time(NULL));
    _Complex double** objVal = *cMatCpy(simDim(2, dim, dim), matrix);
    _Complex double** objVec = cNewM(SHAPE_SQUARE, 'U').create_sm(dim);
    _Complex double** tempMat = cNewM(SHAPE_RECTANGLE, 0).create_rm(simDim(2, dim, dim));
    for (int i = 0; i < dim-2; ++i) {
        _Complex double* tmpHouseV = (_Complex double*)calloc(dim-i-1, sizeof(_Complex double));
        for (int j = 0; j < dim-i-1; ++j) {
            tmpHouseV[j] = objVal[i+j+1][i];
        }
        for (int j = 0; j < dim-i-1; ++j) {
            for (int k = 0; k < dim-i; ++k) {
                tempMat[j][k] = objVal[j+i+1][k+i];
            }
        }
        for (int j = 0; j < dim-i-1; ++j) {
            for (int k = 0; k < dim-i; ++k) {
                double tempSum = 0.;
                for (int l = 0; l < dim-i-1; ++l) {
                    tempSum += tmpHouseV[l]*tempMat[l][k];
                }
                objVal[j+i+1][k+i] -= 2.*tmpHouseV[j]*tempSum/cLenVec(dim-i-1, tmpHouseV)/cLenVec(dim-i-1, tmpHouseV);
            }
        }
        for (int j = 0; j < dim; ++j) {
            for (int k = 0; k < dim-i-1; ++k) {
                tempMat[j][k] = objVal[j][k+i+1];
            }
        }
        for (int j = 0; j < dim; ++j) {
            for (int k = 0; k < dim-i-1; ++k) {
                double tempSum = 0.;
                for (int l = 0; l < dim-i-1; ++l) {
                    tempSum += tempMat[j][l]*tmpHouseV[l];
                }
                objVal[j][k+i+1] -= 2.*tempSum*tmpHouseV[k]/cLenVec(dim-i-1, tmpHouseV)/cLenVec(dim-i-1, tmpHouseV);
            }
        }
        for (int j = 0; j < dim-i-1; ++j) {
            for (int k = 0; k < dim-i; ++k) {
                tempMat[j][k] = objVec[j+i+1][k+i];
            }
        }
        for (int j = 0; j < dim-i-1; ++j) {
            for (int k = 0; k < dim-i; ++k) {
                double tempSum = 0.;
                for (int l = 0; l < dim-i-1; ++l) {
                    tempSum += tmpHouseV[l]*tempMat[l][k];
                }
                objVec[j+i+1][k+i] -= 2.*tmpHouseV[j]*tempSum/cLenVec(dim-i-1, tmpHouseV)/cLenVec(dim-i-1, tmpHouseV);
            }
        }
        freeCVec(tmpHouseV);
    }
    _Complex double delta = 1.e+24 - 1.e+9i;
    _Complex double* retVals = cNewV(0).create(dim);
    do {
        for (int i = 0; i < dim; ++i) {
            retVals[i] = objVal[i][i];
        }
        _Complex double shift = objVal[dim-1][dim-1];
        _Complex double** mh = cNewM(SHAPE_RECTANGLE, 0).create_rm(simDim(2, 3, 3));
        mh[0][0] = objVal[0][0]*objVal[0][0] + objVal[0][1]*objVal[1][0] - 2.*shift*objVal[0][0] + pow(cabs(shift), 2.);
        mh[1][0] = objVal[1][0]*(objVal[0][0] + objVal[1][1] - 2.*shift);
        mh[2][0] = objVal[1][0]*objVal[2][1];
        for (int i = 0; i < 3; ++i) {
            mh[i][0] = 1.*rand()/RAND_MAX;
            mh[i][1] = 1.*rand()/RAND_MAX;
        }
        cSchmidt(3, mh);
        for (int i = 0; i < dim; ++i) {
            tempMat[0][i] = objVal[0][i];
            tempMat[1][i] = objVal[1][i];
            tempMat[2][i] = objVal[2][i];
        }
        for (int i = 0; i < 3; ++i) {
            for (int j = 0; j < 3; ++j) {
                objVal[i][j] = 0.;
                for (int k = 0; k < 3; ++k) {
                    objVal[i][j] += mh[i][k]*tempMat[k][j];
                }
            }
            for (int j = 3; j < dim; ++j) {
                objVal[i][j] = 0.;
                for (int k = 0; k < 3; ++k) {
                    objVal[i][j] += mh[i][k]*tempMat[k][j];
                }
            }
        }
        for (int i = 0; i < dim; ++i) {
            tempMat[i][0] = objVal[i][0];
            tempMat[i][1] = objVal[i][1];
            tempMat[i][2] = objVal[i][2];
        }
        for (int i = 0; i < 3; ++i) {
            for (int j = 0; j < 3; ++j) {
                objVal[i][j] = 0.;
                for (int k = 0; k < 3; ++k) {
                    objVal[i][j] += tempMat[i][k]*mh[j][k];
                }
            }
        }
        for (int i = 3; i < dim; ++i) {
            for (int j = 0; j < 3; ++j) {
                objVal[i][j] = 0.;
                for (int k = 0; k < 3; ++k) {
                    objVal[i][j] += tempMat[i][k]*mh[j][k];
                }
            }
        }
        for (int i = 0; i < dim; ++i) {
            tempMat[i][0] = objVec[i][0];
            tempMat[i][1] = objVec[i][1];
            tempMat[i][2] = objVec[i][2];
        }
        for (int i = 0; i < 3; ++i) {
            for (int j = 0; j < 3; ++j) {
                objVec[i][j] = 0.;
                for (int k = 0; k < 3; ++k) {
                    objVec[i][j] += tempMat[i][k]*mh[j][k];
                }
            }
        }
        for (int i = 3; i < dim; ++i) {
            for (int j = 0; j < 3; ++j) {
                objVec[i][j] = 0.;
                for (int k = 0; k < 3; ++k) {
                    objVec[i][j] += tempMat[i][k]*mh[j][k];
                }
            }
        }
        for (int i = 0; i < dim-3; ++i) {
            double* tmpHouseV = simDVec(3, objVal[i+1][i], objVal[i+2][i], objVal[i+3][i]);
            for (int j = i+1; j < i+4; ++j) {
                for (int k = 0; k < dim; ++k) {
                    tempMat[j][k] = objVal[j][k];
                }
            }
            for (int j = i+1; j < i+4; ++j) {
                for (int k = 0; k < i+1; ++k) {
                    double tempSum = 0.;
                    for (int l = 0; l < 3; ++l) {
                        tempSum += tmpHouseV[l]*tempMat[l+i+1][k];
                    }
                    objVal[j][k] -= 2.*tmpHouseV[j-i-1]*tempSum/dLenVec(3, tmpHouseV)/dLenVec(3, tmpHouseV);
                }
                for (int k = i+1; k < i+4; ++k) {
                    double tempSum = 0.;
                    for (int l = 0; l < 3; ++l) {
                        tempSum += tmpHouseV[l]*tempMat[l+i+1][k];
                    }
                    objVal[j][k] -= 2.*tmpHouseV[j-i-1]*tempSum/dLenVec(3, tmpHouseV)/dLenVec(3, tmpHouseV);
                }
                for (int k = i+4; k < dim; ++k) {
                    double tempSum = 0.;
                    for (int l = 0; l < 3; ++l) {
                        tempSum += tmpHouseV[l]*tempMat[l+i+1][k];
                    }
                    objVal[j][k] -= 2.*tmpHouseV[j-i-1]*tempSum/dLenVec(3, tmpHouseV)/dLenVec(3, tmpHouseV);
                }
            }
            for (int j = 0; j < dim; ++j) {
                for (int k = i+1; k < i+4; ++k) {
                    tempMat[j][k] = objVal[j][k];
                }
            }
            for (int j = 0; j < i+1; ++j) {
                for (int k = i+1; k < i+4; ++k) {
                    double tempSum = 0.;
                    for (int l = 0; l < 3; ++l) {
                        tempSum += objVal[j][l+i+1]*tmpHouseV[l];
                    }
                    objVal[j][k] -= 2.*tempSum*tmpHouseV[j-i-1]/dLenVec(3, tmpHouseV)/dLenVec(3, tmpHouseV);
                }
            }
            for (int j = i+1; j < i+4; ++j) {
                for (int k = i+1; k < i+4; ++k) {
                    double tempSum = 0.;
                    for (int l = 0; l < 3; ++l) {
                        tempSum += objVal[j][l+i+1]*tmpHouseV[l];
                    }
                    objVal[j][k] -= 2.*tempSum*tmpHouseV[j-i-1]/dLenVec(3, tmpHouseV)/dLenVec(3, tmpHouseV);
                }
            }
            for (int j = i+4; j < dim; ++j) {
                for (int k = i+1; k < i+4; ++k) {
                    double tempSum = 0.;
                    for (int l = 0; l < 3; ++l) {
                        tempSum += objVal[j][l+i+1]*tmpHouseV[l];
                    }
                    objVal[j][k] -= 2.*tempSum*tmpHouseV[j-i-1]/dLenVec(3, tmpHouseV)/dLenVec(3, tmpHouseV);
                }
            }
            for (int j = 0; j < dim; ++j) {
                for (int k = i+1; k < i+4; ++k) {
                    tempMat[j][k] = objVec[j][k];
                }
            }
            for (int j = 0; j < i+1; ++j) {
                for (int k = i+1; k < i+4; ++k) {
                    double tempSum = 0.;
                    for (int l = 0; l < 3; ++l) {
                        tempSum += objVec[j][l+i+1]*tmpHouseV[l];
                    }
                    objVec[j][k] -= 2.*tempSum*tmpHouseV[j-i-1]/dLenVec(3, tmpHouseV)/dLenVec(3, tmpHouseV);
                }
            }
            for (int j = i+1; j < i+4; ++j) {
                for (int k = i+1; k < i+4; ++k) {
                    double tempSum = 0.;
                    for (int l = 0; l < 3; ++l) {
                        tempSum += objVec[j][l+i+1]*tmpHouseV[l];
                    }
                    objVec[j][k] -= 2.*tempSum*tmpHouseV[j-i-1]/dLenVec(3, tmpHouseV)/dLenVec(3, tmpHouseV);
                }
            }
            for (int j = i+4; j < dim; ++j) {
                for (int k = i+1; k < i+4; ++k) {
                    double tempSum = 0.;
                    for (int l = 0; l < 3; ++l) {
                        tempSum += objVec[j][l+i+1]*tmpHouseV[l];
                    }
                    objVec[j][k] -= 2.*tempSum*tmpHouseV[j-i-1]/dLenVec(3, tmpHouseV)/dLenVec(3, tmpHouseV);
                }
            }
            freeDVec(tmpHouseV);
        }
        _Complex double*** giv = cCreateGivensMatrix(dim, simDim(2, dim-2, dim-1), cabs(objVal[dim-2][dim-2]), simDVec(2, carg(objVal[dim-2][dim-2]), carg(objVal[dim-1][dim-1])));
        cGivensTransDirect(dim, &objVal, cCvtM().trans(simDim(2, dim, dim), *giv));
        for (int i = 0; i < dim; ++i) {
            for (int j = 0; j < dim; ++j) {
                tempMat[i][j] = objVec[i][j];
            }
        }
        for (int i = 0; i < dim; ++i) {
            for (int j = 0; j < dim; ++j) {
                objVec[i][j] = 0.;
                for (int k = 0; k < dim; ++k) {
                    objVec[i][j] += tempMat[i][k]*(*giv)[j][k];
                }
            }
        }
        freeCMat(dim, *giv);
        freeCMat(3, mh);
        delta = 0.;
        for (int i = 0; i < dim; ++i) {
            delta += cabs(retVals[i] - objVal[i][i]);
        }
    } while (cabs(delta)/dim > eps);
    freeCMat(dim, tempMat);
    return retVals;
}
